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Kinetic  Theories  for  Biofilms 


Qi  Wang  ‘and  Tianyu  Zhang  ' 


Abstract 

We  apply  the  kinetic  theory  formulation  for  binary  complex  fluids  to  develop  a  set  of  hydrodynamic 
models  for  the  two-phase  mixture  of  biofilms  and  solvent  (water).  It  is  aimed  to  model  nonlinear  growth 
and  transport  of  the  biomass  in  the  mixture  and  the  biomass-flow  interaction.  In  the  kinetic  theory 
formulation  of  binary  complex  fluids,  the  biomass  consisting  of  EPS  (Extracellular  Polymeric  Substance) 
polymer  networks  and  bacteria  is  coarse-grained  into  an  effective  fluid  component,  termed  the  effective 
polymer  solution;  while  the  other  component,  termed  the  effective  solvent,  is  made  up  of  the  collective 
ensemble  of  nutrient  substrates  and  the  solvent.  The  mixture  is  modeled  as  an  incompressible  two- 
phase  fluid  in  which  the  presence  of  the  effective  components  are  quantified  by  their  respective  volume 
fractions.  The  kinetic  theory  framework  allows  the  incorporation  of  microscopic  details  of  the  biomass 
and  its  interaction  with  the  coexisting  effective  solvent.  The  relative  motion  of  the  biomass  and  the 
solvent  relative  to  an  average  velocity  is  described  by  binary  mixing  kinetics  along  with  the  intrinsic 
molecular  elasticity  of  the  EPS  network  strand  modeled  as  an  elastic  dumbbell.  This  theory  is  valid  in 
both  the  biofilm  region  which  consists  of  the  mixture  of  the  biomass  and  solvent  and  the  pure  solvent 
region,  making  it  convenient  in  numerical  simulation  of  the  biomass-flow  interaction.  Steady  states  and 
their  stability  are  discussed  under  a  growth  condition.  Nonlinear  solutions  of  the  three  models  developed 
in  this  study  in  simple  shear  are  calculated  and  compared  numerically. 


1  Introduction 

Biofilms,  consisting  of  myriad  microbes,  their  excretions,  and  trapped  particles,  are  ubiquitous  in  nature, 
medical  implants,  rusty  pipes,  and  dentistry  etc.,  where  microbes  survive  on  wet  surfaces.  In  principle, 
a  biofilm  community  can  be  formed  by  a  single  bacterial  species  in  damp  environment;  but,  in  nature, 
biofilms  almost  always  consist  of  rich  mixtures  of  multiple  species  of  bacteria  as  well  as  fungi,  algae, 
yeasts,  protozoa,  other  microorganisms,  debris  and  corrosion  products  etc.  Biofilms  arc  held  together  by  a 
network  of  sugary  molecular  strands  produced  by  the  microbes,  collectively  termed  “extracellular  polymeric 
substances”  or  “EPS”.  In  bacterial  biofilms,  bacterial  cells  arc  held  together  by  the  network  consisting  of 
the  EPS  strands  produced  by  the  bacteria,  allowing  them  to  develop  complex,  three-dimensional,  resilient, 
attached  communities  [10,  11,  13,  28], 

The  center  for  disease  control  and  national  institute  of  health  estimated  that  65%  to  85%  of  all  chronic 
infections  can  be  attributed  to  bacterial  biofilms  [11],  In  human  diseases,  biofilm  infections  are  some  of 
the  most  recalcitrant  to  treat.  Even  with  rigorous  antibiotic  regimens,  some  biofilms,  such  as  those  within 
the  thick  airway  mucus  of  cystic  fibrosis  (CF)  patients,  persist  throughout  the  course  of  the  disease  process 
[18],  One  noted  that  the  gene  expression  of  a  single  bacterial  (planktonic)  cell  differs  drastically  from 
that  of  the  biofilm  colonies  indicating  environmentally  induced  genetic  change  to  the  bacterial  cell.  This 
is  one  of  the  reason  why  antibiotics  that  arc  effective  to  treat  single  bacterium  may  not  be  effective  in  the 
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treatment  of  the  bacteria  encased  in  biofilms.  Bacterial  biofilms  can  also  be  utilized  in  bio-terrorism  in  which 
persistent  ‘bio-terrorist  agent  biofilms'  with  Francisella  tularensis  can  grow  on  surfaces  where  environmental 
amoeba  can  phagocytose  them,  allowing  for  growth  of  fibrosis  [18],  Biofilms  cost  the  U.S.  literally  billions 
of  dollars  every  year  in  energy  losses,  equipment  damage,  product  contamination  and  medical  infections. 
Understanding  the  mechanism  for  the  dynamical  growth,  transport,  detachment  and  break-down  of  biofilms 
is  important  for  improving  water  treatment,  medical  treatment  of  diseases,  protecting  equipment  or  device 
from  corrosion  and  contamination,  and  even  preventing  bio-terrorism.  It  can  have  a  profound  impact  on 
environmental  sciences,  medicine,  civil  engineering,  naval  sciences,  military  applications  and  homeland 
security. 

Modeling  biofilms  has  been  a  challenge  given  their  complex  cellular  structures  and  changing  genetic 
dynamics  in  the  presence  of  foreign  agents  [10,  13,  28].  Biofilms  essentially  behave  like  biogels  in  which 
small  solvent  and  nutrient  molecule  can  permeate  the  network  formed  by  the  bacterial  and  the  EPS  strands.  A 
set  of  discrete  or  semi-discrete  models  associated  with  cellular  automata  have  been  used  to  study  the  various 
aspect  of  the  biofilm  phenomena  in  multidimensional  settings  [29,  30,  31,  32].  Recently,  there  have  emerged 
a  host  of  continuum  models  for  biofilms  treating  them  as  mixtures  of  multiple  idealized  species  [9,  21,  34,  1, 
38,  39].  In  these  models,  the  mixture  is  either  modeled  as  a  multi-fluid  mixture  [9,  21,  34,  1]  or  a  single  fluid 
of  multi-components/species  [38,  39].  In  the  multifluid  models  [9, 21,  1],  which  one  uses  in  the  development 
of  hydrodynamical  models  for  hydrogels,  the  velocity  for  each  individual  component  is  introduced  and  the 
balance  laws  for  the  mass  and  momentum  of  each  species  are  preserved.  The  technical  difficulty  working 
with  these  multi-fluid  models  is  that  the  individual  velocity  fields  of  the  species  arc  used  as  primary  field 
variables  in  these  theories;  however  they  arc  essentially  immeasurable  in  experiments  since  hydrodynamic 
measurements  are  often  limited  to  the  mean  or  average  quantities  in  multiphase  fluids.  Moreover,  the  inflow 
and  outflow  boundary  conditions  for  the  individual  species  are  elusive  when  it  comes  to  mathematically 
solving  the  governing  system  of  equations  without  additional  simplifications.  The  single  fluid  multi-species 
model  however  can  overcome  this  technical  difficulty  by  providing  explicit  relative  velocities  in  terms  of 
the  excessive  velocities  generated  by  nonequilibrium  thermodynamics  (mixing  dynamics)  with  respect  to 
an  average  bulk  velocity  [38,  39]  and  have  been  successfully  used  in  modeling  biofilm  expansion/growth, 
shedding  or  streaming,  detachment  of  biofilm  blobs  from  large  colonies,  and  rippling  phenomenon  under 
shear  [39].  In  this  single  fluid  multi-component  fluid  system,  the  inter-penetration  of  different  species  is 
carried  out  by  their  respect  entropy  and  the  mixing  free  energy.  This  molecular  mechanism  compensates 
the  so-called  friction  forces  due  to  the  relative  motions  among  the  various  species  commonly  seen  in  the 
multi-fluid  theories  near  hydrodynamic  equilibrium.  By  devising  appropriate  intermolecular  potential  for 
the  multiphase  fluid,  we  can  tailor  our  models  to  account  for  the  physics  of  mixing  effectively. 

The  mathematical  models  for  hydrodynamics  of  biofilms  available  so  far  arc  primarily  coarse-grain 
models  in  which  little  cellular  kinetics  is  incorporated.  There  arc  some  common  ingredients  in  the  mul¬ 
tiphase  fluid  models  for  biofilms  (either  multi-fluid  or  single  fluid  ones)  [9,  1,  21,  34,  38],  in  which  the 
nutrient  substrate  is  passively  treated  as  a  part  of  the  solvent  and  the  volume  of  the  bacteria  is  collectively 
treated  a  part  of  the  effective  polymer.  The  effective  polymer  adopted  here  is  also  collectively  identified  as 
the  biomass  in  the  literature.  In  this  paper,  we  interchangeably  use  the  name  of  effective  polymer  or  biomass 
in  order  to  make  contact  with  the  complex  fluid  models  to  be  developed  as  well  as  the  biological  community 
where  biomass  is  the  more  familial-  term. 

The  volume  fraction  of  the  biomass  §n  and  that  of  the  solvent  (f>iS  are  two  basic  hydrodynamic  variables 
along  with  the  velocity  of  the  biomass  v„  and  the  velocity  of  the  solvent  vv  in  multifluid  models  or  the 
average  velocity  v  in  single  fluid  multicomponent  models,  the  concentration  of  the  nutrient  substrate  c,  the 
concentration  of  the  bacteria  B,  the  pressure  /;,  and  the  polymeric  stress  for  the  effective  polymer.  The 
volume  fraction  of  the  biomass  and  the  effective  solvent  are  assumed  to  obey  reaction-transport  equations 
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in  these  models, 


Q+V-Q«V«)=gn, 

%+V-{fas)=g„ 


(1-1) 


where  gn  and  gs  is  the  polymer  production  rate  and  the  solvent  consumption  rate,  respectively.  The  bacterial 
production  as  well  as  the  substrate  consumption  equation  arc  given  as  follows: 


f  +V-(Bvn)=gb, 

I  (4>.vc)  +  V  •  (cv^  -  A4,Vc)  =  -gc, 


(1.2) 


where  gi,  is  the  production  rate  for  the  bacteria,  gc  is  the  consumption  rate  of  the  nutrient  substrate  in  the 
solvent,  Ds  is  the  diffusion  constant  for  the  nutrient. 

The  binary  theories  for  biolilms  differ  in  the  formulation  of  various  velocities  in  conservation  of  mass 
and  momentum  equations  and  the  constitutive  relations  for  gn  and  gs.  Some  end  up  with  a  quasi-compressible 
constraint  for  an  average  velocity  [9,  21,  1]  while  others  retain  the  incompressible  constraint  on  the  average 
velocity  field  [38].  Besides  these  differences,  the  constitutive  equation  for  the  stress  and  reactive  rates  arc 
essentially  the  same  when  both  the  effective  solvent  and  the  effective  polymer  arc  treated  as  viscous  fluids 
which  is  perfectly  valid  in  the  growth  time  scale  of  the  biofilm.  In  this  case,  the  extra  stress  tensor  and  some 
of  the  growth  as  well  as  decay  rates  are  given  by  the  following  constitutive  laws. 


T;  —  2r|„D„,  T ,  —  2]],I), . 

gc  =  pj,  gn  =  ~^+c  ’ 


(1.3) 


where  q,l  v  arc  the  viscosity  of  the  effective  polymer  and  the  solvent,  respectively,  D,1V  =  ^  [Vv„iiS  + 
is  the  rate  of  strain  tensor  for  the  effective  polymer  and  the  solvent,  respectively,  A  is  the  maximum  con¬ 
sumption  rate  of  the  nutrient  substrate,  /./  is  the  maximum  production  rate  of  the  biomass,  Kq  and  Kc  are  two 
half-saturation  constants.  We  note  that  the  bacterial  concentration  decouples  from  the  rest  of  the  equations 
in  the  binary  biofilm  models.  It  can  thus  be  ignored  completely  in  the  models  when  the  focus  is  on  the 
growth  of  the  biomass  or  effective  polymer. 

In  the  phase  field  model  we  developed  in  [38],  the  ensemble  of  the  bacteria  is  effectively  modeled  as 
the  “viscous  solvent”  which  blends  with  the  EPS  network  to  function  effectively  as  a  polymer  solution.  A 
distribution  function  for  the  polymer  network  strands  can  then  be  introduced  to  describe  the  concentration 
of  the  polymer  network  strand  in  the  ’’viscous  solvent”  as  well  as  the  coarse-grain  or  meso-structure  of  the 
polymer  network.  This  thus  motivates  us  to  seek  a  kinetic  formulation  of  biofilm  theories  to  describe  the 
EPS  network  immersed  in  the  viscous  bacterial  bath  and  thereby  to  guide  the  study  of  the  biomass-flow 
interaction. 

In  this  paper,  we  refine  our  single  fluid  multicomponent  formulation  of  biofilm  theories  by  developing 
a  set  of  kinetic  theory  models  for  biofilm- solvent  mixtures  systematically  [2],  This  new  formulation  of  the 
biofilm  models  provides  a  well-structured  theoretical  framework  for  expanding  the  theory  to  incorporate 
additional  molecular  information  and  biochemical  details  of  the  biofilm  components  into  the  model  as  one’s 
understanding  of  biofilm  dynamics  progresses.  In  what  follows,  we  treat  the  biomass-solvent  mixture  as 
incompressible  with  a  divergence-free  average  (bulk)  velocity  field  and  then  develop  a  kinetic  theory  for 
the  effective  polymer  and  solvent  mixture  taking  into  account  the  polymer  conformational  dynamics  of  the 
EPS  strand,  production  of  biomass  at  the  microscopic  level,  and  the  biomass-solvent  interaction.  Within 
this  framework,  a  variety  of  molecular  models  can  be  crafted  to  account  for  the  polymer  network  strand 
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(dumbbell,  FENE  dumbbell.  Rouse  chain,  etc.)  and  a  spectrum  of  polymer  and  solvent  interaction  dynamics 
can  be  proposed  (like  mixing  theory  of  Flory-Huggins,  etc.).  A  simplified  mechanical  model,  the  dumbbell 
model,  for  the  EPS  network  strand  and  a  kinetic  network  theory  along  with  the  Flory-Huggins  mixing 
dynamics  is  set  up  to  illustrate  the  idea  and  demonstrate  how  consistent  stress  expressions  for  the  effective 
polymer  can  be  derived.  The  kinetic  theories  developed  here  arc  de-facto  phase  field  models  with  the  field 
variable  naturally  given  by  the  effective  polymer  (or  biomass)  volume  fraction  (| )„(x,f)  (the  zeroth  moment 
of  the  effective  polymer  distribution  function),  in  which  the  biofilm-solvent  interface  is  given  by  the  zero- 
level  surface  of  the  phase  variable:  lim§_>0+{x|c|)„(x,t)  =  8}  [35,  36].  The  biofilm  region  is  presented  by  the 
non- vanishing  biomass  volume  fraction  variable  (j>„.  In  this  formulation,  we  tacitly  assume  the  EPS  density 
is  directly  proportional  to  the  bacterial  density  and  thereby  a  single  volume  fraction  can  be  employed  to 
represent  the  biomass.  This  phase  field  model  provides  a  low  (computational)  cost  mathematical  formulation 
of  the  complex  interfacial  problem. 

The  rest  of  the  paper  is  organized  as  follows.  First  we  develop  a  set  of  kinetic  theories  for  biofilms  by 
accounting  for  transport  of  polymer  network  strands  and  their  conformational  changes  as  well  as  nutrient 
substrates,  and  the  response  of  the  polymer  network  in  flow  in  three  plausible  ways  within  the  theoretical 
framework  of  network  theories  for  one  fluid  multi-component  mixtures.  We  then  analyze  the  stability  of 
steady  states  to  identify  the  unstable  modes  in  long  wave  regimes.  Finally,  we  study  the  biofilm  expan¬ 
sion/growth  in  one  space  dimension  numerically  and  compare  the  predictions  made  with  the  three  different 
kinetic  theories. 

2  Kinetic  Theories  for  Biohlms 

We  model  the  biofilm  as  a  fluid  mixture  of  two  effective  components:  the  biomass  as  the  effective 
polymer  solution  including  both  the  viscoelastic  EPS  polymer  network  and  the  viscous  bacterial  ’’solution”, 
and  the  effective  solvent  consisting  of  the  solvent  and  the  dissolved  nutrient  substrate.  In  this  theory,  we 
neglect  the  mass  of  the  nutrient  and  instead  only  account  for  its  reactive  effect.  The  EPS  polymer  network 
is  described  by  network  strands  connected  at  junction  points.  The  dynamics  of  the  network  consist  of 
the  dynamics  of  the  center  of  mass  of  the  polymer  strand  and  the  dynamics  for  the  center  of  mass.  This 
description  clearly  cast  a  two-scale  view  on  the  nature  of  the  polymer  network  dynamics.  The  dynamics 
for  the  center  of  mass  of  the  polymer  strand  can  be  described  at  the  macroscopic  or  the  continuum  level; 
whereas  the  dynamics  of  the  polymer  strand  relative  to  the  center  of  mass  belongs  to  the  mesoscale.  Based 
on  this,  we  develop  a  kinetic  theory  for  the  fluid  mixture  accounting  for  the  microstructure  conformation 
of  the  polymer  network  strand  as  well  as  the  interaction  between  the  effective  components,  in  which  we 
model  the  strand  in  the  EPS  network  using  an  elastic  dumbbell  model  and  bacteria  as  viscous  fluids.  In  this 
setting,  the  effective  polymer  is  effectively  treated  as  a  polymeric  solution,  where  the  EPS  polymer  network 
is  “immersed”  in  the  “bacterial  solution.” 

We  denote  the  statistical  weight  of  a  polymer  network  strand  with  its  center  of  mass  located  at  x  and  the 
end-to-end  vector  q  at  time  t  as  \|/(x,  q,  t )  such  that  the  polymer  network  volume  fraction  (|>,2  is  defined  as  the 
zeroth  moment  of  the  statistical  weight  \\i  with  respect  to  the  configuration  variable  q, 

<Mx,0  =  [  V^q=(l),  (2.1) 

Jr3 

where 

((•))=  /,  Wvrfq-  (2-2) 

Jr3 

Notice  that  wherever  the  region  is  filled  with  pure  solvent,  i.e.,  (|)/f (x. r)  =  0,  v|/(x,q,t)  =  0.  Wherever 
(j)„(x.f)  /  0,  however,  a  probability  density  function  can  be  identified  and  a  probability  ensemble 


4 


average  with  respect  to  the  conformation  variable  q  can  be  defined  by 


<(•)>,= 


<(•)>■ 


(2.3) 


To  make  a  distinction  between  the  two  ensemble  averages  defined  above,  we  remark  that  (}q  is  a  probability 
ensemble  average  while  ()  is  an  ensemble  proportional  to  the  volume  fraction  (|).  Apparently,  the  effective 
polymer  volume  fraction  §n (x,  t)  plays  the  role  of  a  phase  field  variable  in  the  theory.  I.e.,  when  §„  (x,  t)  =  0, 
the  fluid  consists  of  entirely  the  solvent;  otherwise,  it  is  a  true  binary  mixture.  Therefore,  the  resulting 
theory  serves  as  an  effective  phase  field  model.  The  two  distinctive  phases  are  differentiated  by  (|)„(x./j  =  0 
(solvent)  and  (| )„(x,f)  >  0  (biofilm  mixture),  respectively.  For  biofilms,  (|)„(x.t)  <  1  is  normally  assumed, 
excluding  the  very  unlikely  situation  where  the  biofilm  is  completely  dry  without  any  solvent. 

In  this  model,  we  treat  the  entire  material  system  as  an  incompressible  single  fluid  mixture  of  two 
effective  components,  in  which  an  average  velocity  v  is  assumed  to  exist  and  divergence  free.  We  propose 
that  the  free  energy  density  of  the  mixture  system  is  composed  of  the  mixing  free  energy  for  the  solvent  and 
the  polymer  network,  elastic  energy  for  the  polymer  network,  self-energy  of  the  polymer  and  the  solvent 
respectively,  and  the  conformational  entropy  for  the  polymer  network  strand  (elastic  dumbbell): 


F 


—  fint  T  fee  +  fes  +  fsol , 


(2.4) 


where  fmt  is  the  extended  Flory-Huggins  mixing  free  energy  density,  fec  is  the  entropic  and  the  elastic 
energy  density  for  the  conformation  of  the  polymer  strand,  fes  is  the  elastic  energy  for  the  polymer  network, 
and  fsoi  is  the  self-energy  for  the  solvent.  Specifically, 

fee  =  vkBT  [  [ln(-|j-)  +^||q||2]t|t(x,q,t)(jfq,  (2.5) 

Jr3  <Pk 

where  fyjksT  is  the  stiffness  of  the  polymer  chain,  v  is  the  number  density  of  the  polymer  chain  and  (In  ^)\|/ 
corresponds  to  the  entropic  contribution  of  the  polymer  chain,  which  is  understood  to  be  possibly  nonzero 
only  in  the  region  where  §n  /  0  and  zero  elsewhere;  the  extended  Flory-Huggins  mixing  free  energy  density 
fint  as  a  function  of  is  given  by 


fint  =kBT 


(2.6) 


where  yi  and  72  measures  the  strength  of  the  entropic  conformational  and  bulk  mixing  free  energy,  respec¬ 
tively,  x  is  the  Flory-Huggin’s  mixing  parameter,  N  is  the  generalized  polymerization  index  for  the  polymer 
strand; 


fes  =  Kf(  B),  (2.7) 

where  fes  is  a  prescribed  function  of  the  Finger  tensor  B  whose  transport  equation  is  given  by 

+  V  •  (v„B)  +  W„  •  B  —  B  •  W „  —  ci\  [D„  •  B  +  B  •  D„]  +  D„  =  2G,Dn,  (2.8) 

where  "kj  is  the  relaxation  time,  —I  <a\  <  I  is  a  model  parameter,  v„  is  the  network  velocity  (defined 
below),  W„  and  D„  are  the  vorticity  and  rate  of  deformation  tensor  associated  to  the  network  deformation, 
respectively,  a\  is  a  parameter  between  -1  and  1,  T|/>  is  a  viscosity  parameter  and  G,  =  ^  is  the  elastic 
modulus; 


fsol  —  (1  §11  )Csol  ‘ 


(2.9) 
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where  Cso\  is  the  self-energy  of  the  solvent  which  is  a  constant.  Since  fsoi  is  a  linear  function  of  the  volume 
fraction,  its  contribution  to  the  chemical  potential  is  merely  a  constant  and  therefore  can  be  neglected.  We 
will  thereby  drop  the  self-energy  of  the  solvent  from  the  free  energy  functional  to  simplify  the  presentation. 
We  note  that  Y2  =  V/V  is  proportional  to  the  reciprocal  of  the  volume  of  the  solvent  molecule  and  ||  •  ||  denotes 
the  I2  norm  of  a  vector.  The  entropic  conformational  free  energy  is  included  in  the  Flory-Huggins  mixing 
free  energy  to  minimize  the  spatial  inhomogeneity  in  the  biomass. 

One  choice  of  the  nonlinear  elastic  energy  for  the  deformation  of  the  center  of  mass  of  the  polymer 
strand  in  the  network  is  given  by 


/(B)=^/r(B)  +  ^lndet(B), 


(2.10) 


where  £,2  and  ^3  are  two  model  parameters.  Then, 


df( «) 

3B 


(2.11) 


Given  the  free  density  functional,  we  can  calculate  the  chemical  potential  of  the  polymer  system  as 
follows 


_  6F  _  S/fe  8fec  fyes 

^  5\| 1  3\\i  3\\i  3\\i  ' 


where 

^=^kBT\\q\\2+vkBT[\nll 

%  =  ^  =  ~kBTylA^n-y2kBT  [^(1 +ln(^)) +ln(l -c^„)  +  1 -X  +  2%^]  , 


S fa 

S\j/ 


/(B). 


(2.12) 


(2.13) 


We  note  that  the  identity:  ^  =  1  is  used  in  the  above  derivation.  So, 


H  =  kBT 


¥ 


-1 


v^ll  q  1 1  +  V  (In  — )  -  Yi  A<j)„  -  y2  [ — (ln(<|)„ ) )  +  In  ( 1  -  <>„ )  +  2%  <!>„] 


N 


+  f(B)  +  const.  (2.14) 


With  the  chemical  potential,  we  next  derive  a  set  of  kinetic  models  for  the  biofilm  fluid.  We  will  present  three 
distinct  versions  by  postulating  three  plausible  ways  that  the  conservative  force  in  translational  diffusion  of 
the  EPS  network  can  be  approximated  in  the  mean-field. 


2.1  Separable  Model  1  for  Biofilms 

In  the  biomass,  we  assume  (1).  the  creation  rate  and  the  annihilation  rate  of  polymer  network  strands  is 
balanced  by  a  growth/production  rate,  (2).  the  translational  diffusion  is  carried  out  through  a  configurational 
space  averaged  chemical  potential,  (3).  the  deformation  of  the  EPS  strand  (elastic  dumbbell)  in  the  config¬ 
urational  space  (q)  is  due  to  the  polymer  velocity  gradient  V v„  to  be  derived  below.  The  transport  equation 
for  the  statistical  weight  \| f  is  then  given  by  the  following  Smoluchowski  equation  [3,  15] 

| V  +  V •  (v¥)  =  V •  (/V (/i)9y)  +  yq- 1 (' V9Ju¥)  -  V  (( Vv*  +  (a -  1  )D„)  •  qy)  +  G« , 

(2.15) 

—  VrKc+c, 
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where  X  is  the  mobility  coefficient,  C,  is  the  friction  coefficient  for  the  rotational  motion  of  the  polymer  chain, 
a  is  a  slip  parameter  between  -1  and  1  accounting  for  slippage  between  the  polymer  and  the  solvent,  Vv„ 
is  the  velocity  gradient  tensor  with  respect  to  the  polymer  velocity,  D„  =  |  [Vv„  +  Vv^]  is  the  rate  of  strain 
tensor  with  respect  to  the  polymer  velocity,  Gn  is  the  polymer  network  strand  production  rate  in  unit  volume 
in  phase  space  (x,q),  /ur  is  the  maximum  growth  rate  of  the  biomass  and  K,  is  a  half-saturation  constant  for 
the  nutrient  consumption  kinetics  via  a  Michaelis-Menton  model.  We  note  that  both  k  and  C,  can  be  second 
order  tensors  when  the  diffusion  is  anisotropic.  In  this  paper,  we  adopt  an  isotropic  diffusion  coefficient 
tensor  for  illustration  purposes  though.  To  simplify  the  presentation,  we  stipulate  that  the  ensemble  average 
((•))<?  =  0  whenever  (|)„  =  0  throughout  the  paper. 

2.1.1  Transport  Equations  For  Volume  Fractions  And  Structure  Tensors:  Low  Order  Moment  Equa¬ 
tions 

Since  polymer  volume  fraction  (])„  is  the  zeroth  moment  of  \\i,  we  obtain  the  transport  equation  of  (|)„  by 
integrating  the  Smoluchowski  equation  over  the  configurational  space  (q-space), 

§L+V.(^v)  =  V-[HlV(q)9]+g„, 


d )nc 

8n=Mr^rc, 


(f*)q 


-kBT  yiA<])„ 


yiksT 


N 


+  ln(l  —  +  2%(|),i 


+ 


(2.16) 


v^kBT(\\q\\2)q+vkBT (ln(^))9  +  /(B)  +  const,  in  biofilm. 

We  remark  that  the  mobility  vanishes  in  pure  solvent  and  so  does  §n{n)q  =  0.  So,  this  is  a  singular  or 
modified  Cahn-Hilliard  equation  [4,  5]. 

From  the  transport  equation  for  the  volume  fraction,  the  local  instantaneous  velocity  can  be  identified 
as  consisting  of  two  parts:  the  average  velocity  v  and  the  excessive  velocity  due  to  the  binary  mixing.  The 
latter  contribution  to  the  flux,  which  can  be  read  off  from  the  transport  equation,  is  assumed  proportional  to 
the  mixing  force  given  by  the  gradient  of  the  ensemble  average  of  the  chemical  potential 

ven  =  -XV(v)q.  (2.17) 

It  is  called  the  excessive  polymer  velocity  that  is  only  defined  in  the  biofilm  region.  The  polymer  velocity  is 
then  identified  as 


V„=v  +  vfr  (2.18) 

This  polymer  velocity  is  not  a  divergence-free  vector  although  the  average  velocity  v  may  be.  In  fact,  this 
velocity  coincides  with  the  average  velocity  outside  the  biofilm  region;  so  it  is  a  globally  defined  hydrody¬ 
namic  variable.  Using  this  notation,  the  transport  equation  for  (j)„  simplifies  to 

^+V'(<|>nV„)=g„.  (2.19) 

With  this  polymer  velocity,  we  can  define  the  rate  of  strain  tensor  and  vorticity  tensor  associated  with  the 
velocity  field: 


D„  =  2(Vv»  +  Vv„),W„  =  —  (Vvn  VVj 


(2.20) 
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We  now  introduce  the  second  moment  of  q  with  respect  to  \\r. 


Q  =  (qq)  =  <Mqq)?-  (2-21) 

The  transport  equation  for  the  second  order  tensor  Q  is  obtained  from  the  Smoluchowski  equation  (2.15)  by 
multiplying  it  with  qq  and  then  integrating  it  over  the  configurational  space  (q): 

|q  +  v  •  (v„Q)  -  W„  Q  +  Q  •  W„  -  a[ D„  Q  +  Q  D„]  =  (<M  -  2^Q)  +  gn Q,  (2.22) 

where  gn  = 

When  (j)n  /  0,  we  introduce  the  second  moment  of  q  with  respect  to  a  probability  density  function 

M  =  (qq)9,  Q  =  <|>„M.  (2.23) 

The  transport  equation  for  M  is  given  by 

(\>n  5-M  +  v„  •  VM  —  W„  •  M  +  M  •  W„  —  a[D„  •  M  +  M  •  D|;]  =  2v^bT  (I  -  2$M).  (2.24) 

at  J  C, 

Canceling  (|);,  /  0  ,  we  arrive  at 

+  v„  •  VM  —  W„  •  M  +  M •  W„  —  a [D„  •  M  +  M •  D„]  =  2V^(I-2^M).  (2.25) 

We  emphasize  that  this  equation  is  valid  only  if  <|)„  /  0.  For  simulation  puiposes,  the  quantity  Q  is  well 
defined  globally  and  easier  to  handle  since  it  vanishes  in  the  pure  solvent  region  along  with  (j>„  =  0. 

We  next  assume  the  statistical  weight  \|/  is  separable,  i.e.,  it  is  a  product  of  the  pdf  \| /„  in  the  configura¬ 
tional  space  of  the  dumbbell  chains  and  the  volume  fraction  ()„: 

\|/(x,q,f)  =(|)n(x,f)\|/„(x,q,f).  (2.26) 

Then,  the  average  chemical  potential  reduces  to 

{n)q  =  -kBT^n-l2kBT  [- ^  +  ln(  1  - <>„)  +  2%^)„]  +vtyBT(\\q\\2)q  +  vkBT(lnyn)q  +  f(B).  (2.27) 

Combining  the  Smoluchowski  equation  and  equation  (2. 19),  we  arrive  at  the  equation  for  the  configurational 
pdf  \| fn  at  4>„  /  0 

=  T  Vn  •  V(V„)  =  v9  •  V^  •  ((Vvn  (//  1  )Dn)  '  qt|tH) .  (2.28) 

It  can  be  readily  verified  that  the  second  moment  equation  of  the  pdf  \| fn  is  given  exactly  by  (2.25). 

Notice  that  the  incompressibility  condition  (|)iV  +  (\)H  =  1  requires 

+V  -  (<M)  =  -XV-  (i,V(q)9)  -g„.  (2.29) 

From  this  transport  equation  for  (j)s,  the  excessive  solvent  velocity  can  be  identified  as 

>l  =  rS), 

Y.v 
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(2.30) 


The  actual  solvent  velocity  is  calculated  by 


W  =  v  +  v*.  (2.31) 

Then,  we  arrive  at  that  the  average  velocity  v  is  in  fact  the  volume  averaged  velocity 

v  =  f,v„  +  ^.v.s.  (2.32) 

Hence,  the  transport  equation  for  (j>v  can  be  recast  into 

^  +  V-(^v,)  =  -g„.  (2.33) 

The  nutrient  substrate  is  assumed  to  be  transported  along  with  the  solvent  velocity.  Its  transport  equation 
is  given  by 


^  (4>vc)  +  V  •  (cv,4>,  -  Adis  Vc)  =  -gc,  (2.34) 

where  c  is  the  nutrient  concentration  per  unit  volume  of  solvent,  dbc  is  the  actual  concentration  per  unit 
biofilm  volume,  and  gc  =  —  is  the  consumption  rate  for  the  nutrient.  Here  A  is  a  maximum  constant 

decay  rate  and  ka  is  the  half-saturation  constant  in  this  model. 

When  the  solvent  is  modeled  as  a  viscous  fluid,  the  constitutive  equation  for  the  extra  stress  in  solvent 
is  given  by 


x.v  =  2r|,D.s,  (2.35) 

where  rp  is  the  solvent  viscosity  and  Dv  =  \  [Vv.s  -P  VvJ ]  is  the  rate  of  strain  tensor  with  respect  to  the  solvent 
velocity. 

Given  the  transport  equations  for  the  volume  fractions,  the  elastic  deformation  tensor  B,  and  the  second 
order  tensor  Q,  we  next  derive  the  elastic  stress  tensor  to  couple  the  meso-structural  variables  (()„.  Q)  and 
the  continuum  variable  B  to  the  macroscopic  momentum  transport  of  the  fluid  mixture. 


2.1.2  Constitutive  Stress  Equations  For  Effective  Polymers 

The  extra  stress  for  the  polymer  in  the  mixture  will  supply  the  crucial  link  to  close  the  governing  system 
of  equations  for  the  biofilm  model.  Given  the  composition  of  the  biomass,  its  contribution  to  the  stress 
tensor  consists  of  two  parts:  the  viscous  stress  §nxps  due  to  the  viscous  bacteria  and  the  elastic  stress  due  to 
the  EPS  network.  The  bacterial  viscous  stress  is  given  by 


xps  =  2rp,iVD„,  (2.36) 

where  v\p  is  the  bacterial  viscosity.  We  use  a  virtual  work  principle  to  calculate  the  elastic  stress  of  the 
polymer.  We  define  the  free  energy  of  the  mixture  system  as 


A  = 


(2.37) 


where  the  free  energy  density  F  is  given  by  (2.4).  The  variation  of  the  free  energy  density  is  calculated  as 
follows 


5  F 


8(/i«r  +  fes)  ,  8/ec  8 fes 

__5+„  +  _8V+w.5b. 


(2.38) 
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We  use 


§<K  =  =  -V  •  (v„<|>n)&, 

5y  =  |\|f8 1  =  -V  •  (v„y)&  -  Vq  •  ((Vv„  +  (a  -  1)D„)  •  q\j/)  &,  (2.39) 


5B  =  §8t  =  [-  [V  •  (y„B)  +  W„  •  B  —  B  •  W„  —  a\  [D„  •  B  +  B  •  D„]]  +  2GjD„]8f , 

assuming  that  the  virtual  time  scale  8t  is  so  small  that  the  virtual  deformation  dominates  the  other  motions. 
It  then  follows  that 


5a  =  /  ^g^(-V  •  (vA))5  tdx  -  ff  [V  •  (v„V)  +  V,  •  ((Vv„  +  (a  -  1)D„)  •  qtfJStdxdq- 


8\)/ 


/  %  :  [[V  •  (v„B)  +  W„  •  B  -  B  •  W„  -  oi  [D„  B  +  B  D„]]  -  2G,Dn]8r<fe. 


(2.40) 


From  the  above  expression,  we  identify  the  elastic  force  as 

F e  =  +  If )  -  (V^)  -  (V|f )  :  B  =  -Yt kBTV  •  (V^V^)  -  4>„V/(B)  -  VPl  -  V^B^lAl) 

where  p\  is  a  scalar  functional  of  (j>„,  and  the  elastic  stress  as 

Gi-Ob  8/es 

-  • - r 

(2.42) 


T  _  G+ll/v  §£ 

Ip  — 


_ / V7  °Jec  n\  i  G~l)  /„T7  &fec\  i  Gl  +  1 )  §fes  1>  i  O Jes  I 

Pt  \ V  <?  5y  9/  2  \(Iv?5y/“h  2  5B  C  +  2  D  '  SB  T 


Combining  Fe  and  xe,  the  extra  elastic  stress  tensor  for  the  effective  polymer  is  expressed  as 

Tf1  =  -Y1^rV^V^+2^vW„M  +  ai^-B  +  G;-(^  +  (^)r).  (2-43) 


Here  we  assume  fes  is  made  up  of  all  bulk  terms.  If  we  introduce 

T  =  2^vkBr(M-^),  (2.44) 

then,  x  satisfies  the  following  transport  equation, 

3 — b  v«  •  Vx  —  Wn  •  T  +  T  •  W„  —  a[D„  •  x  +  T  •  D„]  +  t—  =  D„.  (2.45) 

dt  A?\  A[ 

where  7.|  =  4^kgT  is  the  relaxation  time,  T|„  =  is  the  EPS  polymer  viscosity  [3].  The  rubber-elastic 
model  can  be  viewed  as  a  limiting  case  of  the  current  model  as  A ]  — >  «=  and  a  =  1.  Since  laqvknT (j)„M  = 
M |)„T  +  avknT t[)„I,  by  absorbing  avkuT c^I  into  the  pressure  term,  the  total  extra  stress  is  given  by 

% total  —  t)vxv  T  tytAps  T  Tn  —  t).sX.s'  YikgTV(|)nV(])n  T  cu j)nx  +  tyn^ps  T  nj  ^  •  B 

(2.46) 


After  combining  the  terms  that  can  be  identified  as  a  paid  of  the  stress,  the  elastic  force  is  left  with 

Fe  =  -^V/(B)-V^^:B-Vp1.  (2.47) 

In  summary,  the  kinetic  theory  for  biofilms  consists  of  four  sets  of  equations. 
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Momentum  and  continuity  equation 


V  •  v  =  0, 

Pit  =  ^ V?  +  +  fli  laf  -B  +  Gj(^  +  (^f)r))  (2.48) 


-  V7;  +  y1^rV-(V^Vf1)  +  ^V/(B)  +  V%:B  . 

Here  the  velocity  is  approximated  as  solenoidal  and  the  hydrodynamic  pressure  p  contains  various  scalar 
functionals  of  (])„  and  B  as  well  as  the  static  pressure. 

Transport  equation  for  the  nutrient 

|  ((f>,sc)  +  V  •  (cv,(f>,  -  DsfysVc)  =  -gc.  (2.49) 

Transport  equation  for  the  effective  polymer  volume  fraction 

+  V '  (<!>«▼)  =  v  •  [H.V(/i>9]  +gn.  (2.50) 

Constitutive  equations  for  stress  tensors 

f  +  v„  •  Vx  -  W„  •  x  +  x  •  W„  -  a [D„  •  x  +  x  •  D„]  +  ^  =  ^D„ ,  ips  =  2r|/)D„ ,  x,  =  2r|,D, , 

(2.51) 


f  +  V  •  (v„B)  +  W„  B  -  B  W„  -ai  [D„  B  +  B  D„]  +  £  =  2G,D„. 


2.2  Separable  Model  2  for  Biofilms 


Instead  of  using  the  ensemble  averaged  chemical  potential  to  calculate  the  force  in  the  spatial  transport 
of  \| /,  we  adopt  an  averaged  force  (' Vp)q .  The  Smoluchowski  equation  for  \\i  is  modified  to 

|v  +  V  •  (w)  =  v  •  (A,(V/i)9y)  +  v?  •  J(v^¥)  -  V9((VvB  +  (a  -  1  )D„)  •  q\| /)  +  Gn.  (2.52) 

The  zero  moment  of  the  pdf  is 

^  +  V-(fJv)  =  V-[MVq)^„]+g,i, 


(Vp)q  =  y[-kBTylA^-y2kBT  +ln(l  -<>„)  +  2%^„  +/(B)]=Vg-, 


(2.53) 


where  F  =  fint+  fes.  Note  that  the  effective  polymer  velocity  is  given  by 


=  v-w£ 


(2.54) 


while  the  solvent  velocity  is 


Vi 


v  + A, 


tyn 

(1  -<M 


(2.55) 


By  taking  the  second  moment  of  q  with  respect  to  the  pdf  \|/,  we  arrive  at  equation  (2.22).  The  extra  stress 
tensor  in  this  model  is  identical  to  the  one  derived  in  the  previous  section.  The  difference  lies  in  the  definition 
of  the  polymer  velocity  v„  though. 

If  we  set  \|t  =  (j)„\|/„,  the  transport  equation  for  \| /„  is  given  exactly  by  (2.28).  This  version  of  the  model 
without  the  network  center  of  mass  (the  Finger  tensor)  elasticity  was  studied  in  details  in  [38,  39],  where  it 
was  introduced  phenomenologically.  Now,  we  have  shown  that  the  model  can  be  cast  into  a  low  moment 
projection  of  a  kinetic  model. 
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2.3  Nonseparable  Model  for  Biofilms 

The  above  two  models  arc  developed  using  the  separability  of  the  statistical  weight  \\i  under  two  slightly 
different  mean-field  assumptions  in  the  spatial  or  translational  transport.  We  next  consider  the  third  model. 
Instead  of  the  ensemble  averaged  chemical  potential  or  averaged  force,  we  use  the  same  chemical  potential 
in  both  the  translational  and  configurational  space  transport  in  the  Smoluchowski  equation  for  the  transport 
of  \j /, 

|V  +  V  •  (VV)  =  v '  aVpv)  +  Vq  •  VqjLi\\i )  -  V9((Vv„  +  (a  -  1  )D„)  •  q\|f)  +  Gn.  (2.56) 

The  zero  moment  of  the  pdf  with  respect  to  the  configuration  variable  q  is 
^  +  V  •  (4>„v)  =  V  •  [A,(V/i)j  +  gni 


(Vp)  =  ^n{V/u)q  =  §nV[-kBT  YiAc^„- 


(2.57) 


72 kBT  [-^+ln(l-^)  +  2%^j  +/(B)]  =^V^. 

Notice  that  this  transport  equation  is  identical  to  the  one  derived  using  Model  2.  We  take  the  second  moment 
of  q  with  respect  to  \| i, 

|Q  +  V-(v„Q)-W„-Q  +  Q-W„-n[D„-Q  +  Q-D„] 

=  V  •  (X<(V/i  -  <Vp>,)qq»  +  (4>„I  -  2^Q)  +  gn Q 


6  F 


=  v  •  (^(V%qq»  +  ^(<M  -  2^Q)  +  £„Q 
=  V  •  [Xvksr(VQ  -  Vln^Q)]  +  2vkf(<?nl  -  2^Q)  +  £„  Q. 


(2.58) 


Using  Q  =  (|)„M,  it  follows  that 

3 


|M  +  v„-VM-W„  -M  +  M- W„  — a[D„  -M  +  M-D,,] 
l V  •  [^„vk«TVM]  +  (I  -  2^M) 


(2.59) 


When  (])„  /  0,  it  yields 

|M  +  y„  •  VM  -  W„  •  M  +  M  •  W„  -  a{ D„  -M  +  M-D,,]  =  ^V-  [l^„kBT VM]  +  (I  -  2^M) .  (2.60) 

In  this  model,  the  pdf  \|r  is  not  assumed  separable.  Thus,  we  call  it  a  nonseparable  model.  The  elastic  stress 
constitutive  equations  is  dissipative  and  all  stress  constitutive  equations  are  given  by 

-  W„  •  x  +  x  •  W„  -a[ D„  •  x  +  x  •  D„]  +  ^~  =  ^D„  +  +++ v  ■  , Tps  =  2r|/,D,„  x.s  =  2r|,D,(2.61) 

where  X3  =  ^  is  a  parameter  characterizing  the  stress  diffusive  length  scale.  ^  defines  the  stress 

diffusion  coefficient.  In  practice,  especially,  in  numerical  simulation,  we  use  x„  =  (f)„T.  It’s  zero  when 
4>„  =  0.  The  transport  equation  deduced  from  those  of  x  and  (})„. 


■gf  +  V  •  (v„x„)  —  W„  •  x„  +  x„  •  W„  —  a[ D„  •  x„  +x„  •  D„]  +  ■ 
Here  Bi  is  a  elastic  modulus.  (This  needs  to  be  checked.) 


=  XT V '  (Vx»  “  V/ncj+x,,)  +  £„x„<2.62) 
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3  Closure  Approximation  to  Separable  Model  1  and  Remarks 


In  Separable  Model  1  for  biofilms,  the  transport  equation  for  ())„  is  given  by  (2.16).  This  model  couples 
the  Smoluchowski  equation  for  \(/„  to  the  transport  of  the  volume  fraction  <\)n.  In  order  to  decouple  them,  we 
employ  an  approximation  to  the  averaged  chemical  potential  term  (In  V|/„) c/.  We  assume  the  equilibrium  pdf 
of  the  polymer  strand  obeys  the  Gaussian  distribution  when  the  creation  and  annihilation  rate  cancels  each 
other: 


Yn 


y/2 n  VdetM 


lM  .qq 


(3.1) 


We  replace  the  pdf  \| /„  by  this  equilibrium  pdf  in  (In  \\in)q  =  f  In  \(/„\(/„r/q  to  yield: 


ln\|/„\|/;!dq  =  —  -  ln(detM)  +  const. 


(3.2) 


Then, 


vR,  =  v 


-kBT  YiAf, 


yikBT 


^+ln(l-f,)  +  l-X  +  2xf, 


+  /(B) 


+ 


vkBTV  (£trM-  iln(detM)) . 


(3.3) 


With  this  closure  approximation,  the  transport  equation  for  <()„  only  couples  to  the  transport  equation  for  the 
second  moment  tensor  Q  or  the  elastic  stress  T„.  Consequently,  the  Smoluchowski  equation  is  successfully 
decoupled  from  the  governing  system  of  equations  for  the  mixture. 


Remark  1:  Reformulation  of  Stress  Constitutive  Equations 

Analogous  to  (2.62),  we  can  compute  elastic  stress  tensor  x„  =  (])„T  everywhere  in  the  mixture  through 
its  transport  equation  derived  from  the  second  moment  equation  in  the  models  directly.  In  separable  model 
1  and  2,  it  is  given  by 

^  +  V-  (v„T„)  —  W„  •  T„  +  T„  ■  W„  — fl[D„  -T„  +T„  -D„]  +  D„  +  gn1n-  (3-4) 

The  cost  of  computing  this  set  of  equations  is  comparable  to  that  for  the  stress  tensor  t.  For  comparison 
purposes  with  the  results  in  [38],  we  will  focus  our  analysis  and  numerical  results  on  the  set  of  variables 
<])„, X,  v  and  their  transport  equations  in  this  paper. 


Remark  2:  Model  extension 

In  the  aforementioned  derivation,  the  polymer  network  strand  is  modeled  using  the  simplest  molecular 
model,  the  elastic  dumbbell  or  bead-spring  model.  More  realistic  models  such  as  various  FENE  models, 
network  models  and  other  polymer  models  can  be  incorporated  effortlessly.  If  further  information  on  the 
polymer  network  properties  such  as  its  strand  creation  and  breakage  rate  aside  from  the  polymer  production 
due  to  bacteria  is  known  a  priori,  the  polymer  network  dynamics  can  be  described  in  more  details.  When  any 
nonlinear  molecular  model  is  incorporated  into  the  configurational  and  translational  dynamics,  the  resulting 
theory  is  most  likely  nonlinear  so  that  the  Smoluchowski  equation  is  strongly  coupled  to  the  momentum 
transport  equations.  Therefore,  the  cost  of  solving  the  solutions  of  the  system  can  pose  a  new  challenge  if  a 
decoupled  stress  constitutive  equation  is  not  attainable. 

Next,  we  compare  the  predictions  made  by  the  three  biofilm  models  in  one  space  dimension.  Given  the 
similarity  in  the  viscoelastic  equation  for  the  center  of  mass  dynamics  of  the  polymer  strand  and  that  of  the 
polymer  strand  itself  we  neglect  the  viscoelastic  dynamics  for  the  center  of  mass  in  the  rest  of  the  paper. 
This  is  equivalent  to  using  one  elastic  relaxation  tune  constitutive  equation  for  the  biomass. 
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4  Nondimensionalization  and  Boundary  Conditions 


We  use  a  characteristic  time  scale  to,  length  scale  h,  and  the  nutrient  concentration  scale  co  to  nondimen- 
sionalize  the  variables 


t  = 


t 

to 


ph 2  _  xh2 

Jh,Z^Jb 


c 


c 

Co' 


(4.1) 


where  fo  is  a  characteristic  force  scale.  The  following  dimensionless  parameters  arise 


A  i^o/o  p,  Y ikBT  p  Ya kaTIr  p  _  fotp  p  _  /pfp  p  _  fotg 

—  ^4)^1  —  fn  j  A  2  —  i  ties  —  „  i,iy  tie n  —  „  ,.i ,  tie  ns  — 


A 


fo 


n  sh2  > 


_  f)  Dstp 

n„/i2’’'cPs  ~  ij^p>  Us  h2  > 


®  =  g'P  =  B'»<g  +  *"£).  *  =  A<0,  ft  =^0,  ft  =  ft  =  f ,  A,  =  A3  = 


(4.2) 


where  p0  is  an  averaged  density.  Bi  =  1  is  set  to  define  the  characteristic  force  scale  fo  =  as  the  inertia 

‘6 

force.  We  summarize  the  description  of  the  dimensionless  parameters  below: 

•  Res,  Reps  and  Ren  are  the  Reynolds  number  for  the  solvent,  bacterial  “solvent”  and  EPS  polymer, 
respectively. 

•  A  is  the  dimensionless  mobility. 

•  r i  and  IT  measure  the  strength  of  the  conformational  entropy  and  the  bulk  mixing  free  energy. 

•  Bi  is  a  parameter  for  the  inertia,  which  is  set  to  be  1  in  this  study. 

•  Ds  is  the  dimensionless  diffusion  coefficient  for  the  nutrient. 

•  pr  is  the  dimensionless  biomass  maximum  growth  rate  and  A  is  the  dimensionless  decay  rate  for  the 
nutrient. 

•  KC,K0  arc  the  dimensionless  half  saturation  constants. 

•  Ai  is  the  Deborah  number  for  the  EPS  polymer. 

•  p  is  the  dimensionless  density  of  the  mixture. 

•  A-j  =  -rf-  is  the  effective  stress  diffusion  constant. 

For  simplicity,  we  drop  the  ~ on  the  dimensionless  variables  and  the  parameters.  The  system  of  governing 
equations  of  separable  model  2  in  these  dimensionless  variables  is  listed  below 

V- v  =  0, 


P  =  v  •  0»  {ax  +  ^Ds)  -  [Vp  +  Ti  V  •  ( V<|>„  V<j>„ )] , 

|  Ow)  +  V  •  (cv,<|),.  -  DstysVc)  =  -gc, 
f  +V.(M  =  V-[AiV|]+?„, 


At 


|t  +  v„-V(t)-W„-x  +  t-W„ 


nfD„  -x  +  x-D,, 


4"X  “  Re®"’ 


(4.3) 
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where 


8c=^nK^rc,gn=^nK^rc,  V„=V-A V^,  V,  =  V+^V£. 
The  dimensionless  mixing  free  energy  density  is  now  given  by 

/  =  l|v<MI~+r',2  —  In <)>«  +  (1  —  (])«) ln(l — <[>«)  +  x4>« ( 1  —  ^n)  ) 


Analogously,  the  other  dimensionless  equations  can  be  obtained.  To  save  space,  we  will  only  enumerate 
the  ones  that  differ  from  the  corresponding  equations  in  separable  model  2  listed  above.  For  separable  model 
1,  the  governing  system  of  equations  consists  of  all  equations  in  (4.3)  except  that  the  transport  equation  for 
is  replaced  by  the  following: 

^  +  V-(f,v)  =  V-[A  $„V(n)q]+gn, 

(4.5) 

V(/I>,  =  v  [-HA^  -r2  [-^  +ln(l  -^)„)  +2%^]]  +  4V(tr(x)  -  £lndet(x+  £l)). 

The  governing  system  of  equations  in  the  nonseparable  model  is  given  by  (4.3)  except  that  the  constitutive 
elastic  stress  equation  is  replaced  by  by 

Ai  |x  +  v()-V(x)-W„-x  +  x-W„-fl[D„-x  +  x-D„]  +x  =  ^D„  +  ^V- (<t>„Vx).  (4.6) 

The  boundary  conditions  for  the  moments  are  derived  from  the  boundary  condition  for  the  statistical 
weight  \|r  in  the  Smoluchowski  equation  in  each  model.  For  the  three  models,  they  are  respectively. 


n  •  =  0,  (Separable  Model  1)  , 


n  •  =  0,  (Separable  Model  2)  , 


(4.7) 


n  •  Vqtjr  =  0,  (Nonseparable  Model)  . 

where  n  is  the  unit  external  normal  of  the  boundary.  Taking  the  zero  moment  of  the  above  equations,  we 
arrive  at  the  boundary  condition  for  (j)„  in  the  models,  respectively.  For  the  volume  fraction  (|)„,  an  additional 
flux  condition  has  to  be  imposed  in  all  three  models: 

n  •  V<|>«  =  0.  (4.8) 

The  boundary  conditions  for  the  volume  fraction  §n  in  the  nonseparable  model  are: 

n-V(|)„  =  0,  n  •  (Vq)  =  0.  (4.9) 

The  second  boundary  condition  in  the  above  translates  into 

8F 

^n-V^=0-  (4-10) 

0§n 

Combining  with  the  first  boundary  condition  on  (|)„,  we  arrive  at 

^n-V[r1V>„]=0.  (4.11) 
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We  take  the  liberty  to  set 


n-V>„=0.  (4.12) 

Taking  the  second  moment  over  eq.  (4.7-c),  we  arrive  at  the  boundary  condition  for  the  second  moment 
tensor: 


n  •  (Vpqq)  =  0. 

This  can  be  rewritten  into  the  boundary  condition  in  the  elastic  stress  tensor  x: 


IT  8F  r, 

— jWx  +  V— —  (x  +  — ^T) 

N  §(])«  N 


=  0. 


It  translates  into 


(4.13) 


(4.14) 


(|>„n  •  Vx  =  0. 


(4.15) 


Once  again,  we  make  the  choice  to  set 


n-Vx  =  0.  (4.16) 

The  boundary  condition  for  the  nutrient  is  the  flux-free  boundary  condition  at  the  solid  wall  (n  •  Vc  =  0)  and 
the  Dirichlet  boundary  condition  (c  =  c*)  at  any  boundaries  that  have  access  to  a  nutrient  reservoir. 

The  boundary  conditions  in  separable  model  2  are  identical  except  that  there  arc  no  boundary  conditions 
for  the  stress  tensor  x  since  the  stress  constitutive  equation  in  this  model  is  no  longer  diffusive.  The  boundary 
conditions  for  separable  model  1  arc  similar  to  those  of  model  2  except  that  the  second  condition  for  <)„  is 
now  given  by 


=  n 


V  [-IT A4>„]  +  l-  V(tr(x)  -  ^  lndet(x  +  ^1)) 


=  0. 


(4.17) 


5  Steady  States  in  1-D  and  Their  Linear  Stability 

In  this  section  we  examine  the  solution  of  the  governing  system  of  equations  that  depend  on  one  space 
variable  y  e  I  =  [0, 1].  For  the  models  under  the  zero  elastic  stress  and  zero  nutrient  concentration  condition, 
the  three  models  share  the  same  set  of  nontrivial  steady  states  governed  by  the  singular  steady  state  Cahn- 
Hilliard  equation  subject  to  the  normal  flux  free  boundary  condition.  The  set  of  steady  states  is  given  by 


§n  =<|>0,  c  =  0,  x  =  0,  v  =  0. 


(5.1) 


We  linearized  the  equations  about  the  constant  states  and  then  apply  the  normal  mode  analysis.  In  the 
linearized  systems,  the  normal  modes  of  elastic  stress  components  decouple  from  those  of  (])„  and  velocity 
v.  The  decay  rate  in  separable  model  1  and  2  are  identical  and  given  by 


1 

‘A? 


(5.2) 


while  in  the  nonseparable  model  given  by. 


-^-(1  +A3k~), 


(5.3) 
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where  k  is  the  wave  number  of  the  infinitesimal  perturbation.  The  growth  rate  for  (|)„  is  identical  in  all  three 
models  given  by 

— At>0/r2  (r^2  +  r2/"  (<|>o) ) ,  (5.4) 

where  /  =  —  [j^ln(j>„  +  (1  —  (|)„)ln(l  —  ()„)  +  %(|)„  ( I  —  <))„)]  is  the  bulk  mixing  free  energy  density.  The  modes 
for  the  velocity  and  pressure  are  all  decay,  which  can  be  found  in  [38].  If  /"((to)  is  negative,  a  long-wave 

instability  emerges  for  0  <  \k\  <  \J  —  . 

Another  steady  state  solution  family  is  given  by 


4>n  =  0,  c  =  c0,  x  =  o,  V  =  0, 


(5.5) 


where  co  is  a  constant.  The  linearized  stability  analysis  for  this  set  of  solution  shows  a  growth  given  by  gn 
in  all  models  in  all  wave  numbers.  The  corresponding  growth  rate  is  given  by 


JJrC o_ 

kc  +  cq 


(5.6) 


These  two  ’’growth”  mechanisms  dictate  the  biomass  dynamics  in  the  near  equilibrium  state  and  also  impact 
the  nonlinear  dynamics.  The  detail  about  the  derivation  of  the  linearized  growth  rate  for  the  separable  model 
2  can  be  found  in  [38]. 

One  of  the  characteristics  of  the  biofilm  is  its  time-dependent  dynamical  growth.  In  the  rest  of  the  paper, 
we  will  focus  on  the  comparison  of  transient  solutions  of  the  nonlinear  govern  systems  in  the  three  models. 


6  Numerical  Scheme  for  Transient  Biofilm  Models  in  One  Space  Dimension 


In  this  section  we  discuss  the  numerical  methods  used  for  solving  the  nonlinear  systems  of  partial  dif¬ 
ferential  equations  with  emphasis  on  dynamical  growth  and  transport  of  the  biomass  that  is  homogeneous  in 
(x.z,).  namely,  the  1-D  biofilms  with  space  variable  y  €  I  =  [0, 1],  All  unknowns  arc  function  of  (yj)  only. 
We  denote  v  =  (vx,vy,vz),  \en  =  (0.  r’('.0).  Then,  continuity  equation  along  with  the  boundary  conditions 
implies  =  0.  So,  v„  =  v  +  v*  =  (v„*,  vny,vnz)  =  (vx,vey,vz).  Consequently,  vs  =  v  +  vf  =  (vsx,vsy,vsz)  = 
(vx,vesy,vz).  Using  the  boundary  condition  and  the  momentum  balance  equation,  we  arrive  at 


P 


~ r  1  (<t>n,y)2  +  2 


f  1  tyn  3t\v  (L  r) i'„v 

V  Res  c)v  Reps  c)v 


(6.1) 


We  recall  that  separable  model  1  differs  from  separable  model  2  in  that  its  excessive  flux  for  the  transport 
equation  has  an  extra  term  ^A(|)„V  (tr(x)  —  ^  Indct  (t  +  ^1)).  Likewise,  the  nonscparablc  model  differs 
from  separable  model  2  in  the  transport  equation  for  the  extra  elastic  stress  which  has  an  extra  diffusive  term 
•  (<|)„Vt)  on  the  right  hand  side  of  the  stress  evolutionary  equation  in  the  nonscparablc  model. 

In  the  comparative  study,  we  impose  the  following  boundary  conditions  for  separable  model  2: 


U|>— 0  —  0)  Tr|y=l  —  V shear •  A | hi  —  0, 


V^-nU 


v  —  AV  Mr 

o<|)„ 


•n|a/  =  0, 


(6.2) 


[D()5Vc]  •  n|y=0  =  0,  c\y=i  =  c*, 

where  3 1  arc  the  end  points  of  interval  /.  For  the  nonseparable  model,  additional  boundary  conditions  must 
be  imposed  for  the  stress  tensor 

VxXy  •  n|9/  =  VXyy  ■  h  1 37  =  0,  •  •  •  .  (6.3) 
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For  separable  model  1,  one  of  the  boundary  conditions  for  the  volume  fraction  <|)„  (n  •  V3(|)„  =  0  is  replaced 
by 


n  V 


-r!V2<|>„  +  ^(tr(T)-^lndet(T+^I)) 


=  0. 


(6.4) 


The  numerical  scheme  used  to  study  the  nonlinear  dynamics  of  biolilm  growth  is  a  finite  difference 
scheme  with  uniform  spatial  and  time  step  sizes,  denoted  by  Av  and  At,  respectively.  Assuming  interval 
I  =  [0, 1]  is  divided  into  M  uniform  sub-intervals  with  size  Ay  =  1  /M  at  M  +  1  nodes  yo,yi ,  ■  ■  •  ,yM,  we 
denote  the  values  of  the  numerical  solution  of  ^n,c,vx,Zxy,Xyy  at  ( nAtJAy )  by  ^j^v^x^x^W.  Since 
v  •  n|-j/  =  0,  the  discrete  form  of  the  boundary  conditions  (6.2)  and  (6.3)  yield 


€,1  —  C,-l>  C,2  —  §n- 2;  C,M+ 1  —  §n,M-H  §n,M+ 2  ~  §n,M- 2>  C1  ~  C-V  CM  ~ 

y'J  —  o  v"  —v,  r"  ~  Tn  t"  —  r"  t"  ~  Tn  Tn  —  x" 

vx,0  ~  u’  vx,M  ~  v shear,  —  hcy-l’  ^xy,M+ 1  —  lxy,M- 1’  lyy4  —  lyy,~  1’  V,M+ 1  —  lvv.,W  f 


(6.5) 


For  given  solutions  at  time  step  n  —  I  and  /;.  the  polymer  volume  fraction  at  time  step  n+  1,  (|)JJ 
calculated  by  a  9— method  (0  <  9  <  1): 


Yl+ 1 


IS 


_L/'A«+1 
At 


(C+1  -<K)  +  ev-  [v',+1^;;+1  +a^;;+1v  (riV2^+1 +2r2%^+1); 


=  g;i(c+e,c'I+6)-(i-9)v-  [v”^+Acv(r1v2^+2r2xC);  - 


Ar2v- 


^+C+eVln(l-C+9)) 


After  this,  the  volume  fraction  of  the  solvent  at  time  step  n  +  I  is  obtained  by  <|)"+1  =  1  —  (])"+1  and  the 
substrate  concentration  at  time  step  n+  1,  c',+1,  is  calculated  by 


^(C+1c'I+1 


(j>”cn)  -  9V  •  (Ds(|)"+1  Vc”+1  -  ^,+1v;,+1c',+1 


-gc(C+ V+e)  +  (1  -  9)V  •  (DsrsVcn  -  Cv"c”) . 


(6.7) 


The  spatial  discretization  in  all  semidiscretized  equations  is  done  using  central  differences  to  ensure  the 
second  order  accurate  in  space  and  volume  preserving  for  <()„  when  the  growth  is  turned  off.  Here,  the 
extrapolation  is  accomplished  by  (j)|j+e  =  (1  +  9)(|>”  —  9<)J|  1 ,  c”+e  =  (1  +  9)c"  —  9c'!_1  and  the  nonlinear 
functions  gn,gc  and  some  terms  involving  log-function  are  evaluated  at  these  extrapolated  values.  We  use 
9  =  1/2  in  our  simulations,  thus  the  overall  scheme  is  second  order  in  time  and  space. 

The  average  velocity  components  vx  and  the  stress  components  xxy,xyv.  •  •  •  arc  computed  as  follows.  The 
time  discretization  of  the  equation  for  vx  is  given  by 
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(6.8) 


The  spatial  discretization  is  again  central  difference. 

For  separable  model  1  and  2,  all  six  components  of  the  stress  tensor  satisfy  a  generic  equation  of  the 
form 


dZj  j  3X;  j 

-3r+v>-37=fiAVy" 


(6.9) 


Here  F)y(x,v)  represents  the  (i,j)  component  of  the  stress  tensor  and  it  doesn’t  contain  terms  involving 
partial  derivatives  of  x.  Since  =  0  at  y  =  0.1,  there  arc  no  boundary  conditions  for  the  elastic  stress  tensor 
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x  necessary;  thus,  x  actually  satisfies  an  ODE  ^  =  F(x,  v)  at  y  =  0, 1.  Then  at  the  discrete  level,  we  solve 
Xo ,  x^  by  the  following  Runge-Kutta  method. 

x"+1  =t"  +  ^(K1+2K2  +  2K3  +  K4),  (6.10) 

o 

where 

K\  =  F(x",  Vv^),  F2  =  F(x”  +  yF1,V[V^+2V;;+1]), 

A/  v'!-l-v,!+1 

F3  =  F (x"  +  -k2 ,V[  »  2"  D ,  f4  =  f(x" + wf3 ,  Vvf ') . 

We  solve  M  —  1 ,  by  the  following  upwind  scheme 

W  { []  _  sign(»;j+1/2)]v"+1/2(T;+1  -  x;>  + 

[l+sign(v;j_1/2)]vJ3_1/2(Tj  -  xj_,)}  -t-/?(Tj,VvS).  (6.11) 

For  the  nonseparable  model,  we  have  an  extra  diffusive  term  ^-V  •  (c|>„Vx)  in  the  right  hand  side  of 
(6.9)  and  write  it  as  F(x,  Vv„, <])„).  Using  boundary  condition  (6.5),  we  solve  x'j,  0  <  j  <  M,  by  (6.11)  with 
F(x",Vv")  replaced  by  F(x'j,  VvJJ,( ])").  The  derivatives  are  discretized  again  by  central  differences.  The 
overall  scheme  is  second  order  in  space  and  time. 

7  Model  Comparison  in  1-D  Transient  Flows 

Separable  model  2  in  the  viscous  limit  (A,j  — >•  0)  has  been  studied  in  [38,  39,  7]  in  both  1  and  2  space 
dimension.  Here  we  revisit  the  nonlinear  transient  growth  of  the  biofilm  in  1-D  using  all  three  viscoelastic 
models  with  an  emphasis  on  the  comparison  of  the  model  predictions  under  shear.  We  limit  our  studies  to  the 
growth  and  expansion  of  biofilms  that  are  homogeneous  in  the  (x,z)  plane  alone  using  the  numerical  scheme 
alluded  to  in  the  previous  section.  Table  1  lists  the  ranges  and  values  of  the  dimensional  parameters  used 
in  our  numerical  investigations  which  arc  summarized  from  the  currently  available  literature  [23,  9,  21].  In 
the  numerical  study,  the  initial  profile  for  is  chosen  as  a  step  function  (mimicking  a  localized  distribution 
of  the  biofilm  across  the  shear  cell  with  a  sharp  interface  between  the  biofilm  and  the  solvent),  the  initial 
stress  is  assumed  zero,  and  the  nutrient  concentration  is  initially  saturated  at  the  feeding  level.  We  consider  a 
constant  shear  velocity  vT  =  v shear  =  1  imposed  at  y  =  1  and  the  nutrient  concentration  c*  holds  at  a  constant 
value  at  y  =  1.  The  initial  velocity  profile  is  a  nonlinear  smooth  interpolation  between  the  shearing  speed 
and  the  zero  velocity  at  the  lower  boundai'y  y  =  0. 

The  measured  biofilm  relaxation  time  and  viscosity  vary  with  respect  to  the  wall  stress  exerted  on  the 
rheometric  device  used  in  the  measurement  [37,  23].  The  relaxation  time  can  vary  from  4.57  s  to  169.5  s,  an 
order  of  37  folds  difference  between  the  high  and  the  low  value  in  one  study  [37,  23].  Most  theoretical  and 
computational  studies  on  biofilms  have  been  focused  on  the  biofilm  growth  dynamics  which  is  dominated  by 
the  viscous  effect  in  terms  of  biofilm  hydrodynamics.  Since  the  three  kinetic  models  developed  in  this  paper 
differ  predominantly  in  the  viscoelastic  stress  constitutive  equation,  we  conduct  our  comparative  studies  on 
model  predictions  focusing  on  the  viscoelastic  response.  We  use  two  relaxation  time  values  in  which  the 
larger  one  is  55  folds  of  the  smaller  one. 
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Table  1:  Parameter  Values  Used  In  Simulations 


Symbol 

Parameter 

value 

Unit 

T 

Temperature 

303 

Kelvin 

Yi 

Distortional  energy 

1  x  107 

m  kg  s'2 

Y2 

Mixing  free  energy 

1  x  1017 

m”1  kg  s~2 

1 

Flory-Huggins  parameter 

0.55  and  0.65 

Mobility  parameter 

1  x  10-10  ~  1  x  10~8 

kg_1m3s 

N 

Generalized  polymerization  parameter 

1  x  103 

P 

Max.  Production  rate 

1.4  x  10~3 

kgm~3s-1 

Kc 

Half  saturation  constant  for  polymer  growth 

1  x  10~4 

kgm~3 

K0 

Half  saturation  constant  for  nutrient  decay 

5  x  10~4 

kgm~3 

A 

Max.  Consumption  rate 

0.1 

kgnU3s_1 

Ds 

Substrate  diffusion  coefficient 

2.3  x  10^9 

2  —1 

ms  1 

Viscosity  of  the  EPS  network 

0.02 

kgnu's”1 

0/M 

Viscosity  of  the  bacteria  in  the  polymer  network 

4.3  x  102 

kgnu's-1 

0.V 

Viscosity  of  solvent 

1.002  x  10-3 

kgnu's”1 

pn 

Network  density 

1  x  103 

kgm~3 

P.s 

Network  solvent 

1  x  103 

kgm~3 

co 

Characteristic  substrate  concentration 

1  x  10~3 

kgm~3 

h 

Characteristic  length  scale 

2  x  10~4 

m 

to 

Characteristic  time  scale 

40 

s 

Xi 

Relaxation  time 

50  ~  2  x  103 

s 

a 

Slip  parameter  (a) 

0.92 

M 

Number  of  spacial  sub-intervals 

64  ~  256 

7.1  Highly  elastic  regime:  Ai  >>  1 

We  first  examine  the  models  in  the  “highly  elastic”  regime:  Ai  =  2.5  x  103  >>  1  and  0(1). 

Figure  1  depicts  the  profile  of  the  solutions  of  the  three  models  at  t  =  200  and  the  mobility  parameter 
A  =  2.5  x  10  1 .  We  observe  that  the  biomass  grows  due  to  the  bacterium  and  EPS  production  capability  and 
expands  to  the  solvent  region  by  the  excessive  flux  resulted  from  solvent-polymer  mixing  dynamics.  The 
step  profiles  in  Figure  1(a)  arc  at  /  =  0,  and  the  smooth  curves  arc  at  t  =  200.  The  biomass  growth  is 
fueled  by  the  availability  of  the  nutrient  so  that  the  growth  slows  down  when  the  nutrient  concentration  is 
low  and  eventually  stops  when  it  is  depleted.  We  notice  that  the  time  scales  in  both  growth  and  expansion 
are  comparable  at  this  mobility  value  so  that  the  biomass  in  the  biofilm  region  can  be  transported  swiftly  into 
the  solvent  region  to  fuel  the  expansion.  Hence,  we  don’t  see  the  accumulation  of  biomass  at  the  interface 
leading  to  anomalous  growth  there.  The  maintenance  of  sustainable  level  of  nutrient  at  the  interfacial  region 
is  achieved  by  its  proximity  to  the  nutrient  rich  solvent  region  as  well  as  nutrient  diffusion.  The  differences 
among  results  from  the  three  models  arc  essentially  indistinguishable  in  this  comparison.  This  is  expected 
for  separable  model  2  and  the  nonseparable  model  (labeled  as  model  2  and  3  in  the  figures)  since  they  have 
the  same  transport  equation  for  <()„  and  c.  For  separable  model  (labeled  as  model  1  in  the  figures),  since 
the  dimensionless  mobility  parameter  A  is  at  most  2.5  x  10  7  and  very  small,  the  contribution  from  the 
extra  term  (Ap„V(tr(x)  —  A  |ndct(x  +  I) )  in  the  excessive  flux  is  also  very  small.  Thus  the  plotted 
in  Figure  1(a)  shows  no  visible  difference.  A  more  detailed  comparison  is  given  in  Figure  2  which  depicts 
the  difference  between  computed  by  model  1  and  2  (left  picture)  and  by  model  2  and  3  (right  picture) 
at  t  =  200,  respectively.  We  observe  that  model  2  and  3  give  virtually  identical  result  of  <{>„,  which  differs 
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from  the  result  of  model  1  by  an  order  of  1 0-5.  Therefore,  the  coupling  between  and  stress  tensor  x  in 
model  1  does  make  a  difference.  However,  in  the  current  parameter  value  range,  the  difference  is  so  small 
and  essentially  indistinguishable. 

To  reveal  more  details  about  the  nonlinear  expansion  and  growth  of  biomass  and  dynamics  of  other 
components  in  the  mixture,  we  plot  the  profile  of  the  nutrient  concentration  c  at  t  =  200  in  the  figure  as  well. 
The  initial  distribution  of  the  nutrient  is  assumed  uniform  matching  the  boundary  condition  at  the  feeding 
end.  We  observe  that  the  nutrient  concentration  c  depletes  more  quickly  deep  in  the  region  occupied  by  the 
bulk  biomass  than  near  the  biofilm-solvent  interface. 

The  expansion  direction  of  the  biomass  is  perhaps  best  monitored  by  its  velocity.  The  figure  depicts  the 
nonzero  polymeric  velocity  component  vny  =  vey  in  the  redirection  and  the  polymeric  velocity  vnx  =  vx  in 
the  x-direction.  The  velocity  components  in  both  directions  are  fairly  small  in  magnitudes.  Given  the  large 
bacterial  viscosity  in  the  effective  polymer,  the  motion  of  the  biomass  in  the  biofilm  region  is  very  slow 
compared  to  that  of  the  solvent.  The  growth  of  the  biofilm  in  the  y-direction  is  signatured  by  the  largely 
positive  polymeric  velocity  component  vy  although  it  experiences  a  negative  local  minimum  at  the  interface 
at  t  =  200.  The  negative  transient  velocity  indicates  a  transient  decay  in  the  biofilm  at  t  =  200.  In  general, 
an  overwhelmingly  positive  vy  causes  the  slow  growth  of  the  biofilm's  thickness.  The  equation  for  vx  is 
essentially  a  diffusion  equation  with  the  Dirichlet  boundary  condition  at  y  =  0  and  y  =  1.  Hence,  we  see  an 
approximately  linear  profile  in  vx  in  the  regions  occupied  by  the  biofilm  and  the  solvent,  respectively. 

Since  the  EPS  in  the  biomass  is  modeled  as  a  polymer  network  consisting  of  elastic  dumbbell  strands, 
the  elastic  stress  dynamics  is  worthy  of  a  detailed  interrogation  to  explore  the  mechanical  stress  exerted  on 
the  EPS  polymer  network.  Figure  l(e,f)  depict  the  shear  ((f>„Tx>,)  and  the  normal  elastic  stress  component 
( (f>„xvy)  at  /  =  200,  respectively.  Both  stresses  peak  behind  the  interface  between  the  biofilm  and  the  pure 
solvent  in  the  biofilm  region.  The  peak  of  the  shear  and  normal  stress  correlates  well  with  the  peak  of  the 
velocity  component  vy  in  which  there  exists  a  local  maximum  at  the  interface.  The  results  obtained  from  all 
models  are  identical  numerically. 

At  a  smaller  mobility,  the  transient  biofilm  profile  and  other  hydrodynamic  variables  agree  even  better 
among  the  predictions  by  the  three  models  in  the  parameter  regime  investigated.  Figure  3  depicts  the  solu¬ 
tions  for  all  three  models  at  t  =  200  and  A  =  2.5  x  10  9.  The  biomass  volume  fraction  profile  exhibits  some 
active  growth  and  accumulation  of  biomass  at  the  interface  at  the  smaller  mobility  value.  This  active  growth 
takes  a  toll  on  the  nutrient  concentration  at  the  interface  as  well,  where  the  nutrient  concentration  exhibits  a 
local  minimum  corresponding  to  the  local  maximum  in  the  biomass  fraction.  This  is  the  consequence  of  the 
fact  that  the  growth  time  scale  exceeds  that  of  the  transport  time  scale  at  this  mobility  parameter  leading  to  a 
local  accumulation  of  the  biomass  and  thereby  the  noticeable  consumption  of  the  nutrient.  The  elastic  shear 
stress  value  is  considerably  small  in  the  biofilm  than  in  the  previous  case.  The  fluctuation  at  the  interface 
in  the  elastic  shear  stress  is  not  seen  at  t  =  200.  Figure  4  exploits  the  fine  details  of  the  velocity  inside  the 
biofilm  region  during  the  shearing  process  at  /  =  200  and  A  =  2.5  x  10  9.  It  reveals  the  magnitude  of  the 
slow  motion  of  the  biofilm  in  the  x-direction  leading  up  to  the  edge  of  the  interface. 

We  also  examine  the  solution  of  the  system  at  a  different  mixing  parameter  %  =  0.65  when  there  is  a 
growth  mode  due  to  the  polymer-solvent  mixing  dynamics.  Figure  5  depicts  the  solution  at  t  =  200  and 
A  =  2.5  x  10  9 .  When  the  growth  mode  due  to  the  mixing  dynamics  exists,  the  biomass  growth  dominates. 
At  t  =  200,  there  is  little  expansion  of  the  biofilm  into  the  pure  solvent  region  shown.  To  the  contrary,  the 
active  growth  of  the  biomass  tends  to  pull  the  biomass  into  the  the  active  growth  region  around  the  interface 
in  the  cell.  This  is  best  shown  in  the  velocity  vy  where  it  is  negative  at  the  interface.  This  is  related  to  the 
nucleation  process  inherent  to  the  Cahn-Hilliard  dynamics  in  this  parameter  regime.  Meanwhile,  the  elastic 
shear  and  the  normal  stress  component  all  peak  at  the  fastest  growing  spot.  The  situation  is  alleviated  as 
the  mobility  parameter  or  the  shear  speed  increases.  Once  again,  the  three  models  give  virtually  identical 
predictions. 
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(e)  <|>/|Xxy  (f)  4>n Tyy 

Figure  1:  Solutions  of  the  three  models  at  t  =  200,  mobility  parameter  A  =  2.5  x  10-7,  mixing  parameter 
X  —  0.55,  and  relaxation  time  Ai  =  2.5  x  103.  The  separable  model  1  and  2  are  referred  to  as  model  1  and  2 
while  the  nonseparable  model  is  referred  to  as  model  3  in  the  text  and  the  figures,  (a).  The  biomass  volume 
fraction  <j)„;  (b).  the  nutrient  concentration  c;  (c).  the  velocity  component  in  the  x-direction  v*;  (d).  the 
velocity  component  in  the  y-direction  =  v^;  (e).  the  elastic  shear  stress  ;  (f).  the  elastic  normal 
stress  <j)„xvv.  The  three  models’  predictions  are  essentially  indistinguishable  in  this  parameter  regime. 
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Figure  2:  Difference  of  the  biomass  volume  fraction  between  model  1  and  2  as  well  as  model  2  and  3  at 
t  =  200.  (a).  The  difference  in  the  biomass  volume  fraction  between  model  1  and  model  2.  The  largest 
difference  is  in  the  order  of  0(1O-5).  (b).  The  difference  in  the  biomass  volume  fraction  between  model  2 
and  model  3.  Model  2  and  3  give  virtually  identical  results.  Parameter  values  are  identical  to  the  ones  used 
in  Figure  1. 


7.2  Viscoelastic  regime 

We  next  consider  the  regime  where  Ai  decreases  to  50  at  fixed  A  =  2.5  x  10-7.  The  velocity  component 
starts  to  show  deviations  between  the  one  calculated  from  model  1  and  those  by  model  2  and  3  at  t  =  60 
shown  in  Figure  6.  vy  predicted  by  separable  model  2  and  the  nonseparable  model  tends  to  show  tamed 
fluctuations  at  the  interface  than  that  of  separable  model  1 .  The  shear  stress  profile  predicted  by  three  models 
show  some  differences  as  well,  in  that  separable  model  2  and  the  nonseparable  model  give  qualitatively  the 
same  predictions  while  separable  model  1  yields  a  fluctuating  elastic  shear  stress  near  the  interface.  The 
normal  stress  prediction  by  the  three  models  in  <|>nXyy  are  all  comparable  however.  Figure  7  depicts  the 
difference  between  <(>„  computed  using  model  1  and  2  (left  picture)  and  using  model  2  and  3  (right  picture)  at 
t  =  60,  respectively.  We  again  observe  that  model  2  and  3  give  virtually  identical  result  for  <))„,  which  differs 
from  the  result  of  model  1  in  the  order  of  10-3.  Although  the  difference  in  <|)„  is  still  invisible  from  Figure 
6(a),  it  is  two  orders  of  magnitude  larger  than  the  difference  shown  in  Figure  2.  The  apparent  fluctuations 
in  velocity  component  vv  and  thereby  the  elastic  shear  stress  (j)„tvv  from  model  1  are  attributed  to  the  larger 
deviation  in  <(>„  from  model  1  in  the  current  parameter  value  range.  We  remark  that  the  results  shown  in  the 
figure  are  computed  using  mesh  size  Ay  =  1/64,  and  computations  with  mesh  size  Ay  =  1/128, 1  /256, 1/512 
give  the  same  result  which  indicates  that  the  stress  fluctuation  is  space  is  not  an  artifact  of  the  coarse  mesh, 
rather  it’s  the  constitutive  response  to  the  slight  variation  in  the  biomass  volume  fraction  and  the  associated 
excessive  velocity  variation.  Furthermore,  since  the  governing  equations  for  the  stress  x  are  the  same  for 
separable  model  1  and  2,  and  are  solved  by  the  same  numerical  scheme,  the  difference  in  x  should  be 
attributed  to  the  difference  between  these  two  models.  Namely,  <(>„  and  X  are  coupled  in  separable  model  1 
but  not  in  separable  model  2.  We  examine  the  difference  between  separable  model  2  and  the  nonseparable 
model  further  along  in  the  simulation.  At  t  =  200,  the  difference  in  the  elastic  stress  starts  to  emerge  as 
shown  in  Figure  8.  The  solutions  are  comparable  to  those  depicted  in  Figure  1  except  that  the  two  elastic 
stress  components  predicted  using  separable  model  2  differ  from  those  predicted  by  the  nonseparable  model 
in  the  interfacial  region.  The  difference  is  quantitative  though. 

At  mixing  parameter  %  =  0.65,  the  biofilm  dynamics  in  the  viscoelastic  limit  is  analogous  to  the  highly 
elastic  limit;  the  growth  instability  due  to  the  biomass-solvent  mixing  along  with  the  biomass  production 
dominates  leading  to  biomass  growth  in  volume  fraction  in  the  originally  occupied  location  in  biofilm.  There 
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Figure  3:  Solutions  of  the  three  models  at  t  =  200,  mobility  parameters  A  =  2.5  x  10-9,  mixing  parameter 
X  =  0.55,  and  relaxation  time  Ai  =  2.5  x  103.  The  nutrient  is  consumed  the  most  at  the  front  of  the  biofilm. 
The  three  models’  predictions  are  essentially  indistinguishable. 
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Figure  4:  Blow-up  of  v*  in  the  biofilm  region  with  respect  to  solutions  of  all  three  models  at  t  =  200.  The 
results  are  essentially  indistinguishable. 

is  a  slight  pulling  back  at  the  biofilm-solvent  interface  due  to  biomass  nucleation.  The  growth  velocity  is  an 
order  of  magnitude  weaker  than  in  the  highly  elastic  limit.  Figure  9  plots  a  solution  atr  =  200  and  X  =  0.65. 

Finally,  we  comment  on  the  solutions  at  a  smaller  mobility  in  the  viscoelastic  limit.  The  general  trend 
of  the  solution  follows  that  alluded  to  in  the  highly  elastic  limit  for  <(>,,.  The  solutions  predicted  by  the  three 
models  agree  very  well.  The  elastic  shear  stress  are  all  positive  with  the  largest  value  given  by  separable 
model  1. 

In  summary,  the  three  models  give  qualitatively  the  same  results  in  the  parameter  range  we  investigated. 
When  mobility  is  low,  the  growth  mode  dominates.  Otherwise,  the  spatial  expansion  of  the  biofilm  and 
local  growth  occur  simultaneously.  Minor  differences  are  identified  among  separable  model  1 ,  separable 
model  2  and  the  nonseparable  model.  The  volume  fraction  transport  equation  in  separable  model  1  couples 
to  the  stress  constitutive  equation.  Whereas,  that  of  the  separable  model  2’s  equation  for  <)>„  decouples 
from  the  elastic  stress  equation  and  is  driven  by  the  equation  of  c,  so  does  the  nonseparable  model.  But, 
the  nonseparable  model  has  a  diffusive  stress  constitutive  equation  which  couples  to  the  equation  of  <])„ 
and  requires  a  no-flux  boundary  condition  for  the  stress  tensor  at  the  boundary.  Ideally  the  nonseparable 
model  is  the  model  we  should  use  since  it  does  not  impose  any  a  priori  assumption  on  the  distribution  of 
the  EPS  polymer  network  (i.e.,  the  separability  of  the  pdf  distribution).  However,  this  comparative  study 
has  demonstrated  the  competency  of  separable  model  2  in  the  parameter  regime  we  investigated,  which  is 
simpler  than  the  nonseparable  model. 

8  Conclusions 

We  have  systematically  developed  a  set  of  kinetic  theories  for  the  biofilm,  a  mixture  of  biomass  and 
solvent,  using  the  one-fluid  multi-component  formulation  to  model  the  nonlinear  growth  and  transport  of  the 
biomass  (extracellular  polymeric  substances  EPS  and  bacteria)  and  the  interaction  between  the  biomass  and 
nutrient  and  solvent  in  flows.  The  theoretical  framework  allows  detailed  conformational  information  of  the 
EPS  polymer  network  strand  to  be  accounted  for  and  has  the  potential  to  be  expanded  to  incorporate  more 
microscopic  details  about  the  biomass  and  cell-to-cell  communication  like  quorum  sensing  in  the  future. 
Adopting  three  distinctive  formulations  of  the  dominating  mean  field  force  or  velocity  in  the  translational 
diffusion  in  the  Smoluchowski  equation  for  a  normalized  distribution  \|r  (a  statistical  weight,)  we  derive  three 
distinct  models  for  the  biomass-solvent  mixture.  All  these  models  are  valid  not  only  within  the  biofilm  region 
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Figure  5:  Solutions  of  the  three  models  at  t  =  200,  mobility  parameter  A  =  2.5  x  10  9,  mixing  parameter 
%  =  0.65,  and  relaxation  time  Aj  =  2.5  x  103.  The  results  are  essentially  indistinguishable. 
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Figure  6:  Solutions  of  all  three  models  at  t  =  60,  mobility  parameter  A  =  2.5  x  10-7,  mixing  parameter 
X  =  0.55,  and  relaxation  time  Aj  =  50.  Model  1  exhibits  the  largest  fluctuation  in  velocity  and  elastic  stress 
at  the  biofilm-solvent  interface  while  the  other  two  models  show  mild  fluctuations. 
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Figure  7:  Biomass  volume  fraction  difference  between  two  pairs  of  models  at  t  =  60.  (a).  The  difference 
between  model  1  and  model  2  at  t  =  60.  The  largest  difference  is  in  the  order  of  0(  10-3).  (b).  The  difference 
between  model  2  and  model  3  at  t  =  60.  The  two  models  give  virtually  identical  results.  Parameter  values 
are  identical  to  those  used  in  Figure  6. 


Figure  8:  Solutions  of  model  2  and  3  at  t  =  200,  mobility  parameter  A  =  2.5  x  10“7,  mixing  parameter 
%  =  0.55,  and  relaxation  time  A\  =  50.  The  fluctuation  in  stress  is  amplified  in  long  time  demonstrating  that 
model  3  exhibits  the  smallest  fluctuation  among  the  three  models  investigated. 
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Figure  9:  Solutions  of  the  three  models  at  t  —  200  and  mobility  parameter  A  =  2.5  x  10-9.  mixing  parameter 
X  =  0.65,  and  relaxation  time  Ai  =  50.  The  results  from  model  2  and  3  are  essentially  indistinguishable. 
The  stress  fluctuation  from  model  1  differs  slightly  from  those  from  model  2  and  model  3. 
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but  also  in  the  pure  solvent  region  making  them  bona  tide  multiphase  hydrodynamic  models.  A  specific 
contact  is  made  between  separable  model  2  and  the  continuum  phase  field  model  developed  previously 
based  on  a  phenomenological  approach  establishing  the  microscopic  foundation  for  the  phase  held  model, 
which  has  been  used  to  study  biofilm  dynamics  in  3-D  flow  chambers  [38,  7],  We  then  analyze  the  linear 
stability  properties  of  a  set  of  constant  steady  states  shared  by  all  three  models  revealing  the  potential  long¬ 
wave  instability  in  the  models  for  the  biofilm  growth  in  addition  to  the  inherent  biomass  growth  mechanism 
for  all  waves.  Finally,  we  compare  the  three  models  on  their  transient  nonlinear  dynamics  in  1-D  shear  flows 
in  the  viscoelastic  regime  with  only  one  relaxation  time.  Our  numerical  results  show  that  all  three  models 
predict  qualitatively  the  same  behavior  with  the  diffusive  stress  model  rendering  a  smoother  stress  profiles 
and  reduced  stress  values  in  the  parameter  range  investigated.  For  model  1,  the  biomass  profile  and  the 
associated  excessive  velocity  for  the  biomass  leads  to  a  elastic  stress  fluctuation  within  the  biomass-solvent 
interface.  This  feature  is  robust  in  model  1  and  is  believed  to  be  the  result  of  enhanced  coupling  between 
the  stress  and  the  biomass  transport  that  is  unique  in  this  model.  In  the  viscous  limit,  the  three  models 
collapse  into  a  single  multiphase  phase  field  model  which  was  studied  extensively  in  [38,  39,  7],  Extension 
of  this  formulation  can  be  made  to  deal  with  multiple  species  of  biomass  and  thereby  allow  a  full  coupling 
of  quorum  sensing  mechanism.  This  work  will  be  forthcoming  in  a  sequel. 
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